/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
** ****************************************************************** */

//
// 文件: GPUConcreteKernels.cu
// 说明: GPU混凝土本构kernel实现（仅包含kernel函数和C接口）
// 注意: 设备函数实现在GPUConcreteKernels_impl.cuh中
//

// ============================================================
// Performance optimization: Disable kernel debug output
// ============================================================
#ifndef GPU_ENABLE_DEBUG_OUTPUT
#define GPU_ENABLE_DEBUG_OUTPUT 1  // Temporarily enabled for debugging
#endif

#if GPU_ENABLE_DEBUG_OUTPUT
#define KERNEL_DEBUG_PRINT(...) printf(__VA_ARGS__)
#else
#define KERNEL_DEBUG_PRINT(...) ((void)0)
#endif

#include "GPUConcreteKernels.h"
#include "GPUElementSupport.h"
#include <cuda_runtime.h>
#include <stdio.h>

//==============================================================================
// CUDA Kernel实现
//==============================================================================

extern "C" {

/**
 * @brief 初始化材料状态变量kernel
 *
 * 根据单元材料参数初始化每个高斯点的ConcreteState结构体
 * 主要初始化损伤阈值：
 * - rp = ft / E (拉伸损伤阈值)
 * - rn = fc / E (压缩损伤阈值)
 */
__global__ void initMaterialStatesKernel(
    int numGaussPoints,
    const GPUElementData* d_elementData,
    ConcreteState* d_states,
    int numElements,
    int gaussPointsPerElement
)
{
    int gpIdx = blockIdx.x * blockDim.x + threadIdx.x;
    if (gpIdx >= numGaussPoints) return;

    // 确定高斯点所属的单元
    int elemIdx = gpIdx / gaussPointsPerElement;
    if (elemIdx >= numElements) return;

    // 获取单元数据
    const GPUElementData& elem = d_elementData[elemIdx];
    ConcreteState& state = d_states[gpIdx];

    // 初始化损伤变量为0
    state.dp = 0.0;
    state.dn = 0.0;

    // 从单元数据中提取材料参数
    // materialProps布局（根据GPUMaterialExtractor）：
    // [0] = E (弹性模量)
    // [1] = nu (泊松比)
    // [2] = ft (抗拉强度)
    // [3] = fc (抗压强度)
    double E = elem.materialProps[0];
    double ft = elem.materialProps[2];
    double fc = elem.materialProps[3];

    // ✅ 初始化损伤阈值（与本构kernel中的rp0/rn0公式一致）
    // rp0 = ft / sqrt(E)
    // rn0 = sqrt((-k + sqrt(2)) * fc / sqrt(3))
    double f2c = 1.16 * fc;
    double k = sqrt(2.0) * (f2c - fc) / (2.0 * f2c - fc);

    state.rp = (E > 1e-10) ? (ft / sqrt(E)) : 0.0;
    state.rn = (E > 1e-10) ? sqrt((-k + sqrt(2.0)) * fc / sqrt(3.0)) : 0.0;

    // 初始化塑性应变为0
    for (int i = 0; i < 6; i++) {
        state.eps_p[i] = 0.0;
    }

    // 初始化有效应力为0
    for (int i = 0; i < 6; i++) {
        state.sige[i] = 0.0;
    }

    // 初始化辅助变量为0
    for (int i = 0; i < 8; i++) {
        state.stateVars[i] = 0.0;
    }
}

/**
 * @brief 启动材料状态初始化kernel
 */
int initializeMaterialStates(
    int numGaussPoints,
    const void* d_elementData,
    ConcreteState* d_states,
    int numElements,
    int gaussPointsPerElement,
    cudaStream_t stream
)
{
    KERNEL_DEBUG_PRINT("[GPU DEBUG] initializeMaterialStates called: numGaussPoints=%d, numElements=%d, gaussPointsPerElement=%d\n",
                       numGaussPoints, numElements, gaussPointsPerElement);

    if (numGaussPoints <= 0 || d_elementData == nullptr || d_states == nullptr) {
        KERNEL_DEBUG_PRINT("[GPU ERROR] initializeMaterialStates: Invalid parameters (numGaussPoints=%d, d_elementData=%p, d_states=%p)\n",
                           numGaussPoints, d_elementData, d_states);
        return -1;
    }

    int blockSize = 256;
    int gridSize = (numGaussPoints + blockSize - 1) / blockSize;

    KERNEL_DEBUG_PRINT("[GPU DEBUG] Launching kernel with gridSize=%d, blockSize=%d\n", gridSize, blockSize);

    initMaterialStatesKernel<<<gridSize, blockSize, 0, stream>>>(
        numGaussPoints,
        static_cast<const GPUElementData*>(d_elementData),
        d_states,
        numElements,
        gaussPointsPerElement
    );

    // Check for launch errors
    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        KERNEL_DEBUG_PRINT("[GPU ERROR] Kernel launch failed: %s\n", cudaGetErrorString(err));
        return -1;
    }

    // Wait for kernel to complete and check for execution errors
    err = cudaDeviceSynchronize();
    if (err != cudaSuccess) {
        KERNEL_DEBUG_PRINT("[GPU ERROR] Kernel execution failed: %s\n", cudaGetErrorString(err));
        return -1;
    }

    KERNEL_DEBUG_PRINT("[GPU DEBUG] initializeMaterialStates completed successfully\n");
    return 0;
}

/**
 * @brief 测试kernel：批量计算混凝土本构
 */
__global__ void concreteKernel(
    int numGaussPoints,
    const double* d_strainIncr,
    const double* d_strainTotal,
    const PlasticDamageConcreteParams* d_params,
    ConcreteState* d_states,
    double* d_stress,
    double* d_tangent
)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= numGaussPoints) return;

    // 获取当前积分点的数据
    const double* strainIncr = &d_strainIncr[idx * 6];
    const double* strainTotal = &d_strainTotal[idx * 6];
    const PlasticDamageConcreteParams& params = d_params[idx];
    ConcreteState& state = d_states[idx];
    double* stress = &d_stress[idx * 6];
    double* tangent = &d_tangent[idx * 36];

    // 调用本构计算（设备函数在头文件中inline定义）
    computePlasticDamageConcrete3D(
        strainIncr,
        strainTotal,
        params,
        state,
        stress,
        tangent
    );
}

/**
 * @brief 启动混凝土材料kernel
 */
int launchConcreteKernel(
    int numGaussPoints,
    const double* d_strainIncr,
    const double* d_strainTotal,
    const PlasticDamageConcreteParams* d_params,
    ConcreteState* d_states,
    double* d_stress,
    double* d_tangent,
    cudaStream_t stream
)
{
    if (numGaussPoints <= 0) return -1;

    int blockSize = 256;
    int gridSize = (numGaussPoints + blockSize - 1) / blockSize;

    concreteKernel<<<gridSize, blockSize, 0, stream>>>(
        numGaussPoints,
        d_strainIncr,
        d_strainTotal,
        d_params,
        d_states,
        d_stress,
        d_tangent
    );

    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        KERNEL_DEBUG_PRINT("[GPU ERROR] launchConcreteKernel failed: %s\n", cudaGetErrorString(err));
        return -1;
    }

    return 0;
}

} // extern "C"
